SCBI Data

# From forestecology package
# devtools::install_github("rvalavi/blockCV")
library(blockCV)
library(sfheaders)

scbi_2013 <- read.csv("data/scbi.stem2.csv") %>% 
  select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
  mutate(date = lubridate::mdy(date)) %>%
  filter(gx < 300, gy > 300, gy < 600)

scbi_2018 <- read.csv("data/scbi.stem3.csv") %>% 
  select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
  mutate(date = lubridate::mdy(date),
         dbh = as.numeric(dbh)) %>%
  filter(gx < 300, gy > 300, gy < 600)

census_df1 <- scbi_2013 
census_df2 <- scbi_2018
id <- "stemID"

cluster_names <- data.frame(group = c(1:7), group_name = c("high_sla", "low_sla", "shrubs", "tall_light_wood", "evergreens", "oaks", "other"), group_alias = c("High SLA trees", "Low SLA trees", "Shrubs", "Tall trees - Light wood", "Evergreens", "Oaks", "Other"))

spptable <- read_csv("data/scbi.spptable.csv") %>% 
  left_join(cluster_names, by = "group") %>% 
  mutate(group = ifelse(is.na(group), 7, group),
         group_name = ifelse(group == 7, "other", group_name),
         group_alias = ifelse(group == 7, "Other", group_alias))

scbi_growth_df <- 
  # Merge both censuses and compute growth:
  compute_growth(census_df1, census_df2, id) %>%
  # Convert growth from cm to mm to make result comparable
  mutate(growth = growth/10,
             sp = as.factor(sp)) %>% 
  inner_join(spptable, by = "sp") %>% 
  filter(status == "A")
# Add spatial information
cv_fold_size <- 100
max_dist <- 7.5

scbi_study_region <- 
 #tibble(x = c(0,400,400,0,0), y = c(0,0,640,640,0)) %>% 
  tibble(x = c(0,300,300,0,0), y = c(300,300,600,600,3)) %>% 
  sfc_polygon()

# Add buffer variable to data frame
scbi_growth_df <- scbi_growth_df %>%
  add_buffer_variable(direction = "in", size = max_dist, region = scbi_study_region)

scbi_cv_grid <- spatialBlock(
  speciesData = scbi_growth_df, theRange = 100, yOffset = 0.9999, k = 9, verbose = FALSE)

# Add foldID to data
scbi_growth_df <- scbi_growth_df %>% 
  mutate(
    foldID = scbi_cv_grid$foldID
  )

# Visualize grid
scbi_cv_grid$plots +
  geom_sf(data = scbi_growth_df, aes(col=factor(foldID)), size = 0.1)

scbi_cv_grid_sf <- scbi_cv_grid$blocks %>%
  st_as_sf()

Create focal vs. competitor data

focal_vs_comp_scbi <- scbi_growth_df %>% 
  create_focal_vs_comp(max_dist, cv_grid_sf = scbi_cv_grid_sf, id = "stemID")

focal_vs_comp_scbi <- focal_vs_comp_scbi %>% 
  left_join(select(spptable, sp, focal_group = group_name), by = c("focal_sp" = "sp")) %>% 
  left_join(select(spptable, sp, comp_group = group_name), by = c("comp_sp" = "sp"))

Spread into wide form for modeling (note: should we do this by genus?)

focal_vs_comp_scbi_wide <- focal_vs_comp_scbi %>% 
  group_by(focal_ID, focal_group, dbh, foldID, growth, comp_group) %>%
  summarize(comp_basal_area = sum(comp_basal_area)) %>% 
  pivot_wider(names_from = comp_group, values_from = comp_basal_area, values_fill = 0)

group_order <- focal_vs_comp_scbi_wide %>% 
  ungroup() %>% 
  filter(focal_group != "other") %>% 
  count(focal_group, sort = T)

# set factor order by largest group
focal_vs_comp_scbi_wide <- focal_vs_comp_scbi_wide %>% 
  mutate(focal_group = factor(focal_group, levels = c(as.character(group_order$focal_group), "other")))
head(focal_vs_comp_scbi_wide)
## # A tibble: 6 x 12
## # Groups:   focal_ID, focal_group, dbh, foldID, growth [6]
##   focal_ID focal_group   dbh foldID growth high_sla low_sla shrubs
##      <int> <fct>       <dbl>  <int>  <dbl>    <dbl>   <dbl>  <dbl>
## 1        4 low_sla     136.       6  0.103     2.34   4.58    1.96
## 2        5 shrubs       88        6  0.150     1.47   9.88    6.89
## 3       79 tall_light… 477.       9 -0.161    13.8   34.8     4.52
## 4       80 shrubs       51.5      6  0.253     0      0.356   2.45
## 5       96 other        23        9  0.262     5.09   0       1.53
## 6      101 tall_light… 654.       7  0.552     2.83   0.859   0   
## # … with 4 more variables: tall_light_wood <dbl>, other <dbl>, oaks <dbl>,
## #   evergreens <dbl>
# write_csv(focal_vs_comp_scbi_wide, "data/focal_vs_comp_scbi_wide.csv")

Model

model <- lm(growth ~ dbh + focal_group + high_sla + low_sla + tall_light_wood + evergreens + oaks + shrubs + other, data = focal_vs_comp_scbi_wide)

summary(model)
## 
## Call:
## lm(formula = growth ~ dbh + focal_group + high_sla + low_sla + 
##     tall_light_wood + evergreens + oaks + shrubs + other, data = focal_vs_comp_scbi_wide)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57235 -0.06850 -0.01825  0.04713  2.34239 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            8.146e-02  6.274e-03  12.984  < 2e-16 ***
## dbh                    5.346e-04  1.295e-05  41.292  < 2e-16 ***
## focal_grouphigh_sla   -2.108e-02  5.922e-03  -3.559 0.000375 ***
## focal_grouplow_sla    -3.323e-02  6.127e-03  -5.424 6.05e-08 ***
## focal_groupshrubs     -2.605e-03  6.974e-03  -0.374 0.708766    
## focal_groupoaks       -6.129e-02  8.991e-03  -6.816 1.02e-11 ***
## focal_groupevergreens  1.106e-01  1.802e-02   6.134 9.07e-10 ***
## focal_groupother      -5.637e-03  6.095e-03  -0.925 0.355133    
## high_sla              -1.247e-05  1.466e-04  -0.085 0.932230    
## low_sla               -6.453e-05  1.060e-04  -0.609 0.542575    
## tall_light_wood       -3.430e-04  7.381e-05  -4.646 3.45e-06 ***
## evergreens            -8.056e-04  4.491e-04  -1.794 0.072878 .  
## oaks                  -5.625e-04  8.943e-05  -6.290 3.39e-10 ***
## shrubs                 7.497e-03  1.447e-03   5.180 2.29e-07 ***
## other                 -7.219e-04  4.913e-04  -1.469 0.141836    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1413 on 6302 degrees of freedom
## Multiple R-squared:  0.3035, Adjusted R-squared:  0.3019 
## F-statistic: 196.1 on 14 and 6302 DF,  p-value: < 2.2e-16

\[ y = \beta_0 + \beta_1d_f + \beta_2g_f + \sum_{g=1}^G \beta_{g}m_g + \varepsilon \\ \text{ where } d_f \text{ is the diameter at breast height of the focal tree, } g_f \text{ is the cluster group of the focal tree, } \\ g \in \{\text{tall trees/light wood, high SLA trees, low SLA trees, shrubs, oaks, evergreens, other}\}, \\ \text{ and } m_g \text{ is the biomass of all competitor trees of the specific focal group. Coefficients are with respect to tall trees/light wood} \]

Plot all trees in SCBI plot:

scbi_growth_subset <- scbi_growth_df %>% 
  filter(stemID %in% focal_vs_comp_scbi_wide$focal_ID)

ggplot(data = scbi_growth_subset) +
  geom_sf(aes(col = sp), size = 0.5) +
  labs(title = "All trees in SCBI site")

ggplot(data = scbi_growth_subset) +
  geom_sf(aes(col = Genus), size = 0.5) +
  labs(title = "All trees in SCBI site")

ggplot(data = scbi_growth_subset) +
  geom_sf(aes(col = Family), size = 0.5) +
  labs(title = "All trees in SCBI site")

ggplot(data = scbi_growth_subset) +
  geom_sf(aes(col = factor(group_alias)), size = 0.5) +
  labs(title = "All trees in SCBI site")

Plot all the trees, specific species

ggplot(data = scbi_growth_df %>% filter(sp == "nysy")) +
  geom_sf(aes(col = sp), size = 0.5)+
  labs(title = "All Nyssa sylvatica trees in Michigan Big Woods site")

ggplot(data = scbi_growth_df %>% filter(sp == "astr")) +
  geom_sf(aes(col = sp), size = 0.5)+
  labs(title = "All Asimina triloba trees in Michigan Big Woods site")

maples <- read_csv("https://rudeboybert.github.io/SDS390/static/data/maples.csv")

Model loos like this: \[ \begin{eqnarray*} y &=& \beta_0 + \beta_{\text{dbh}}\cdot \text{dbh} + \beta_{\text{evergreen_bm}}\cdot \text{evergreen_bm} + \beta_{\text{maple_bm}}\cdot \text{maple_bm} + \beta_{\text{misc_bm}}\cdot \text{misc_bm}\\ && + \beta_{\text{oak}}\cdot \text{oak} + \beta_{\text{short_tree}}\cdot \text{short_tree} + \beta_{\text{shrub}}\cdot \text{shrub} + \epsilon \end{eqnarray*} \]

where

  1. \(\text{dbh}\) is the “starting” DBH in 2008
  2. \(\text{bm}\) is “biomass” as estimated by the basal area.